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Abstract 

Multifractal analysis studies signals, functions, images or fields via the fluctuations of their local regularity 
along time or space, which capture crucial features of their tenrporal/spatial dynamics. It has become 
a standard signal and image processing tool and is commonly used in numerous applications of different 
natures. In its common formulation, it relies on the Holder exponent as a measure of local regularity, 
which is by nature restricted to positive values and can hence be used for locally bounded functions only. 
In this contribution, it is proposed to replace the Holder exponent with a collection of novel exponents 
for measuring local regularity, the p-exponents. One of the major virtues of p-exponents is that they can 
potentially take negative values. The corresponding wavelet-based multiscale quantities, the p-leaders, are 
constructed and shown to permit the definition of a new multifractal formalism, yielding an accurate practical 
estimation of the multifractal properties of real-world data. Moreover, theoretical and practical connections 
to and comparisons against another multifractal formalism, referred to as multifractal detrended fluctuation 
analysis, are achieved. The performance of the proposed p-leader multifractal formalism is studied and 
compared to previous formalisms using synthetic multifractal signals and images, illustrating its theoretical 
and practical benefits. The present contribution is complemented by a companion article studying in depth 
the theoretical properties of p-exponents and the rich classification of local singularities it permits. 

Keywords: Multifractal analysis, p-exponent, wavelet p-leaders, negative regularity, multifractal 
detrendend fluctuation analysis 


1. Introduction 

Context: Scale invariance and multifractal analysis. The paradigm of scale-free dynamics, scale 
invariance, or scaling, has been shown to be relevant to model empirical data from numerous real-world 
applications of very different natures, ranging from Neurosciences [1, 2, 3, 4], heart rate variability [5, 6, 7], 
bones textures [8] to physics (turbulence [9]), geophysics (rainfalls [10], wind [11], earthquakes [12]), finance 
[13, 14, 15], Internet traffic [16], art investigations [17], music [18] ... In essence, scale invariance is associated 
to signals whose temporal dynamics involve a wide range (continuum) of time/space scales, rather than being 
dominated by one or a few specific time/space scales. 
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Multifractal analysis provides a generic framework for studying scaling properties both theoretically and 
practically. Multifractal analysis describes the dynamics of the fluctuations along time or space of the local 
regularity, h(x o), of the signal (or image) X around location Xq. However, though theoretically grounded on 
local or pointwise quantities, multifractal analysis actually does not aim to produce estimates of h. for all x. 
Instead, it provides practitioners with a global, geometrical and statistical description of the fluctuations of 
h across all time or space locations. This global information is encoded in the so-called multifractal spectrum 
T>(h) (also termed singularity spectrum) which is defined as the fractal dimension of the set of points x 
where the pointwise regularity exponent h(x) takes the value h. 

For real-world data, the estimation of the multifractal spectrum does not rely directly on local estimates 
of h but instead involves a thermodynamics-inspired formalism, the so-called multifractal formalism. A 
multifractal formalism first requires the choice of multiresolution quantities , hereafter labelled T x (a : fc), i.e., 
quantities that depend jointly on the position k and the analysis scale a. Scale invariance is defined as the 
power-law behavior, with respect to the analysis scale a, of the time average of the q —th power of these 
T x (a, k): 

1 n a 

S(a,q) = - y £\ T x(. a ,k)\ q -a«“ ) . (1) 

The quantity «S(a, q) is referred to as the structure function associated with the quantities T\(a , fc), and the 
function f(q) is termed the scaling function. The scaling function f(q) can be related to functional space 
signal or image characterization, cf. e.g., [19, 20]. Most importantly, the multifractal formalism relates the 
Legendre transform C(h) of £(g) to T>(h) via the formula practical estimation of T>(h) can be obtained via 
the: 

V(h ) < C(h) := inf (d + qh — £(<?)) (2) 

(cf. e.g., [19, 21]), and strongly ties the paradigm of scale invariance to multifractal properties and the 
spectrum V(h) and provides a way of estimating V(h) in practice. 

For a number of synthetic processes with known T>(h ), it can be shown that the inequality (2) turns out to 
be an equality. It is also well-known that earlier formulations of multifractal formalisms relying on increments 
or wavelet coefficients as multiresolution quantities yield poor or incorrect estimates of the multifractal 
spectrum. Commonly used formalisms are relying on Wavelet Transform Modulus Maxima (WTMM) [22] 
or MultiFractal Detrended Fluctuation Analysis (MFDFA) (cf. [23] for the founding article, and [24, 25] for 
related developments). Recently, it has been shown that a formalism based on Wavelet Leaders provides 
a theoretically well grounded and practically robust and efficient framework for the estimation of D{h) 
[26, 27, 28, 21], 


Multiresolution quantities and regularity exponents. An often overlooked issue lies in the fact 
that the choice of multiresolution quantity T x is crucial, both theoretically and practically: The multifractal 
formalism (2) relies on the assumption that the pointwise exponent h x can be theoretically recovered through 
local log-log regressions of the multiresolution quantities 


h(x) = liminf 

a —>-0 


log (Tx{a,k(x))) 
log a 


( 3 ) 


(where T x (a, k(x)) denotes the value taken by T x (a , k) at the location k(x) closest to x). Therefore, choosing 
the multiresolution quantity T x (a,k) on which the formalism is based also amounts to choosing (and thus 
changing) the definition of local regularity. The predominantly if not exclusively employed notion of local 
regularity h is the Holder exponent (cf. Section 3.1). Typical examples for T x that have been associated with 
Holder exponents are increments, oscillations, continuous and discrete wavelet coefficients, wavelet transform 
modulus maxima (WTMM) [22, 29, 30, 31], or wavelet leaders [19, 32], Note, however, that oscillations and 
wavelet leaders are the only multiresolution quantities for which it has been shown theoretically that (3) 
holds and that it actually characterizes the Holder exponent [19]. It has been shown that increments and 
wavelet coefficients are associated to another (and weaker) notion of regularity [33, 32, 28], thus partially 
explaining their poor estimation performance. So far, there are no theoretical results available connecting 
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WWTM or MFDFA to Holder regularity. It is shown here in Section 4 below that the latter, MFDFA, is 
not associated to Holder regularity but to a different local regularity measure. At present, wavelet leaders 
constitute the core of the current state-of-the-art multifractal formalism associated to Holder regularity 
[26, 27], 

However, the Holder exponent based measure of regularity suffers from the fundamental drawback that it 
is, by definition, positive and can hence not characterize negative regularity. Multifractal analysis is therefore 
restricted to functions whose local regularity is positive everywhere. This represents a severe limitation since 
data are often found to contain negative regularity in applications (cf. e.g., [28, 20]). 

Goals, contributions and outline. In this context, the present contribution aims to overcome this 
restriction for the application of multifractal analysis by making use of new pointwise regularity exponents, 
the p-exponents, and the corresponding multiresolution quantities, the p-leaders, which have been proposed 
in the mathematical literature recently [34]. The p-exponents offer a more versatile framework than Holder 
exponents and notably enable to theoretically define and measure negative regularity. In the companion 
paper [35], p-exponent regularity and p-leaders have been theoretically defined, studied and characterized. 
The first goal of the present contribution is the development of the corresponding multifractal formalism. 
The second goal is to compare the p-leader multifractal formalism theoretically and practically against 
wavelet leaders and against MFDFA. 

After having formulated a wavelet-based definition of p-exponents and p-leaders and having briefly re¬ 
stated their key properties (Section 2), the p-leader multifractal formalism is derived in Section 3, and 
explicit estimation procedures for the multifractal spectrum that take into account the discrete and finite- 
size nature of data are developed. Section 4 starts with a brief description of the original formulation of 
DFA and recalls its multifractal extension, the MFDFA method. It establishes, for the first time, clear 
theoretical and practical connections between MFDFA and p-exponents and points out the theoretical and 
practical consequences of the conceptual differences between MFDFA and p-leaders. Section 5 provides a 
numerical study and validation of the proposed p-leader multifractal formalism and compares it with the 
wavelet leader and MFDFA formalisms in terms of estimation performance, illustrating the clear practical 
benefits of the proposed novel multifractal formalism. 

A Matlab toolbox implementing estimation procedures for performing p-multifractal analysis will be 
made publicly available at the time of publication. 

2. Wavelet based definition of pointwise p-exponents 

2.1. Wavelet characterization of uniform regularity 

2.1.1. Discrete wavelet decomposition 

Wavelet analysis has already been massively used for the study of the regularity of a function [36, 29, 37, 
38, 19]. Let ... 2 d -i denote a family of mother wavelets, consisting of oscillating functions with 

fast decay and good joint time-frequency localization. Mother wavelets are defined so that the collection of 
dilated (to scales 2 J ) and translated (to space/time location 2 J fc) templates 

{2 for i = 1, - - -2 rf — 1, je Z and k £ Z d }, (4) 

forms an orthonormal basis of L 2 (M. d ) [39]. Without loss of generality, a d-dimensional orthonormal wavelet 
basis is defined here by tensor product of the univariate orthonormal wavelet basis. Mother wavelets 
are moreover defined to guarantee a number of vanishing moments, a strictly positive integer > 1 such 
that J m P(x)if^\x)dx = 0 for any polynomial P of degree strictly smaller than N^, while this integral 
does not vanish for some polynomials of degree A, /; . The coefficients of the d —dimensional discrete wavelet 
transform of X are defined as 

cfl= [ X{x)2- di ^ i) {2~ : >x-k)dx (5) 

J R d 

where the L 1 normalization for the wavelet coefficients is better suited in the context of local regularity 
analysis, see [26, 21]. 
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2.1.2. Wavelet structure function and scaling function 
Let q > 0. The wavelet structure function is defined as 

s c (j, q) = 2* l c Sl 9 > ( 6 ) 

k i —1 


and the wavelet scaling function as 


\/q > 0, 


r](q) = lirninf 

j ->-oo 


log (S c (j, q)) 

log(2J) 


( 7 ) 


The wavelet scaling function iq(q) does not depend on the (sufficiently smooth) wavelet basis which is used, 
see [33, 19] and Section 3.1 of companion article [35] . 

The wavelet scaling function is of central importance for validating uniform regularity assumptions on 
the function X in (13) below, see the companion paper [35] and [20, 33, 21, 28] for details. The wavelet 
scaling function also serves as an important auxiliary quantity in the construction of estimators for finite 
resolution data, cf. Section 3.4. 


2.2. p-leader based definitions of p-exponents 
2.2.1. p-leaders 

Let us now define the multiresolution quantities at the heart of the present contribution, the p-leaders 
(cf. the companion paper [35] for details). Let k = (fci,..., kf), j £ Z and let the dyadic cubes be indexed 
by 

A (= A jik ) := [Vku&ik! + 1)) x ••• x [Vk d ,2?(k d + 1)) 

and, accordingly: c] 1 = c^\. Furthermore, let 3A denote the cube with the same center as that of the cube 
A but 3 times wider. 

Definition 1. Let p > 0 . The p-leaders are defined as [40, 20, 34, 35] 

&\j,k) = e^ = 

where j' < j is the scale associated with 

This definition is illustrated in [35, Fig. 

(13) below. 


2“ —1 


H J2\ c y\ P 2 ~ dU ~ j/) 

A'C3A i =1 


p > o 


( 8 ) 


the sub-cube X' of width 2 J included in 3A. 

1]. The sum in (8) is finite if the condition r fip) > 0 is satisfied, see 


2.2.2. p-exponents 

When p-leaders are well defined (p(p) > 0), the pointwise p-exponent h p (x o) can be defined as follows 
[40, 20, 34], 

Definition 2. Let p > 0 and p(p) > 0. The pointwise p-exponent of X at Xq is 


h p (xo) = lirninf 


j->-oo 


log(2f) 


where Aj(a:o) denotes the dyadic cube of width 2° that contains Xq. 


( 9 ) 
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This definition does not depend on the wavelet as long as the wavelet is smooth enough (cf. [35] Section 3.1 
and [39, Chapter 3]). 

In practice, (9) essentially means that 

4 P) (xo) ~ C2 3h ^ xa) when a = 2 j 0, (10) 

which meets the natural prerequisite (cf. (3)) to construct a multifractal formalism. This will be the subject 
of Section 3. 

When p > 1, Definition 2 is equivalent to the original definition of p-exponents supplied by the T£(xo) 
regularity of A. Calderon and A. Zygmund [41], which essentially states that 

T^\a,x 0 )=[- d ( \X(u + x 0 )-P Xo (u + x 0 )\ p du) ~ Ca h ^ (11) 

\a d J B (0,a) J 


where B(xo,a) denotes the ball of radius r centered at xo, P XQ is a polynomial of degree less than h p (xo) 
and C, R are positive constants, see [34, 40, 20] and [35, Definition 1] for details. Note that the advantage 
of Definition 2 is that it is valid for all p > 0. Furthermore, note that for the case p = +oo, we recover from 
(8) the definition of classical wavelet leaders , see, e.g., [19, 32, 26, 27] 

^A = 4 +00) = Slip l c vl 

iG[l,...,2 d -l], A'C3A 


and from (11) the definition of the classical Holder exponent [19] 


\X(a + Xg) — P x0 ( a + z 0 )| ~ C\a\ h ( x °\ |a| —> 0. 


( 12 ) 


When 0 < p < 1, this wavelet-based Definition 2 of p-exponents has been shown to be well-grounded, 
relevant and useful (cf. [35]). 

2.3. Properties of p-exponents 

The following theoretical key properties are satisfied by p-exponents (details and proofs are given in the 
companion paper [35]). 

Domain of definition. p-exponents are defined for p £ (0,po], where the critical Lebesgue index po is 
defined as 1 

Po = sup{p : r?(p) > 0}. (13) 

In contradistinction, Holder exponents are only defined for locally bounded functions, i.e., functions for 
which po = Too. 

Negative regularity. p-exponents can take values down to —d/p, see [35, Theorem 1], Therefore, p- 
exponents enable the use of negative regularity exponents. For instance, they enable to model local behaviors 
of the form |X(x)| ~ l/|ir — a; 0 | a for a < d/p. 

p-exponents for different values of p. p-exponents for different values of p do not in general coincide 
but satisfy h p (xo) > h p /{x o) for p <p' < po, see [35, Theorem 1]. Therefore, p-exponents for different values 
of p can provide complementary information for the characterization of singularities. 

Classification of singularities. p-exponents enable a generic classification of singularities, see [35, 
Section 4] for details and examples. Let X have a singularity at Xq that has p-exponent h p {X q). The 
singularity of X at Xq is p-invariant if the function p —> h p (xo) does not depend on p. Moreover, 


1 Note that (13) actually provides a lower bound for the quantity po (.xq ) whose definition is given in [35, Theorem 1] and 
will be used here as an operational surrogate. 
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- X has a canonical singularity at £0 if the p-exponent of its fractional integral of order t equals 

h p (xo) + t (cf., [35, Definition 3-4]). Examples are provided by deterministic self-similar functions, e.g., 
the cusps \x — xo|“, a > 0, see [35, Proposition 1], An important property of canonical singularities is 
that they form a subclass of the p-invariant singularities [35, Theorem 2]. 

- A singularity that is not canonical is called an oscillating singularity: 

- An oscillating singularity that is p-invariant is termed balanced , see [35, Definition 7]. An example 
is provided by the “chirp” type singularities \x — xo|“ sin(l/|a: — £o|' 3 ), a,/3 > 0. 

- An oscillating singularity for which the function p —» h p (x 0 ) does depend on p is termed lacunary. 
An example is given by the lacunary comb F ar/ (x) in [35, Eq. (12)] and by the random processes 
in Section 5.1.2. 

Holder exponent. The Holder exponent is given by the p-exponent with p = + 00 : h( xq) = h+^Xo). 
However, it is central to realize that p-exponents for finite values of p measure a regularity in X that does 
not in general coincide with that captured by the Holder exponent, cf., [35, Theorem 1], 

When p —»• 00 , Condition r)(p) > 0 in (13) is the counterpart to the condition met in wavelet leader 
multifractal analysis (cf. [35, 28, 20]): 


log 


h min = liminf 

j -*—00 


sup 

i,k 


AO | 

'j,k\ 


log (21) 


> 0 . 


(14) 


3. Multifractal Analysis 

In practice, the values taken by a pointwise regularity exponent h(x) cannot be extracted point by point 
from data X(a;). Instead, multifractal analysis aims to provide a global description of the fluctuations along 
time or space of h (x), termed the singularity or multifractal spectrum, which can be estimated from data 
X ( x ) by recourse to a so-called multifractal formalism. In the case of Holder regularity, such a formula is 
given by the wavelet leader multifractal formalism, cf., e.g., [19, 32, 26, 27, 21]. The goal of this section 
is to define and study a multifractal formalism based on p-exponents and p-leaders and to derive explicit 
expressions for estimators that can be applied to data in applications. 

The p-leaders and p-exponents have not yet been used for practical multifractal analysis intended for 
real-world data analysis except in our preliminary works [42, 43, 44]. 

3.1. Multifractal p-spectrum 

The multifractal p-spectrum of X is defined as follows. 

Definition 3. The multifractal p-spectrum r D^ p l(h) is the Hausdorff dimension dim// of the set of points 
where the p-exponent takes the value h 

V^{h) = dim// ({a; G M. d , h p (x) = h}) . (15) 

The support of the spectrum is the image of the mapping x —> h p {x), i.e. the collection of values of h such 
that {x £ : h p (x) = h} 7^ 0. 

By convention dim(0) = — 00 . 

Because p-exponents are necessarily larger than —d/p, the support of 'D’A is included in \—d/p, + 00 ]. In 
addition, is further constrained by the following sharper condition: 

Proposition 1. Let p > 0 and let X be a function for which rj(p ) > 0. Then 

Wh< 0, V (p \h) < d + hp. (16) 
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The proof is given in Appendix A. 

This condition implies that the multifractal p-spectrum is necessarily below a straight line that connects 
points (—d/p, 0) and (0, d). This is illustrated in Figs. 1 and 5. 

For further details on multifractal analysis, notably for the definition of the Hausdorff dimension, inter¬ 
ested readers are referred to [45, 37, 38, 19]. 

3.2. p-leader multifractal formalism 

p-leader structure functions. The p-leader multifractal formalism relies on the p-leader structure 
function sj> p \j,q), which is defined as the sample moments of the p-leaders 

VgeK, S ( f\j 1 q) = 2^Y J {^l) q - (17) 

k 

p-leader scaling function. The p-leader scaling function is defined as 

, log (s {p) (j,q)\ 

VgeR, C p (?) = liminf - ^ .... J • (18) 

J^-oo l0g(2J ) 

The function (q) does not depend on the wavelet basis when the mother wavelets are C°° (e.g., in the 
Schwarz class), cf. [34, 40, 19]. 

p-leader multifractal formalism. A multifractal formalism can be defined 
of the mapping q — > (^(q) 

C (p) (h) := inf (d + g/i-C (p) (<?))■ 

It provides practitioners with an upper bound for If r](p) > 0, then 

Wi, f>W(/i)<£W(/i) (20) 

as a particular occurrence of the general principles that lead to (2), see [34, 40, 33]. 

A key property of (20) is that it yields a nontrivial upper bound also for the decreasing part of the 
spectrum (a property that did not hold for formulas based on wavelet coefficients instead of wavelet leaders 
even in the setting supplied by the Holder exponent [19]). This upper bound turns out to be sharp for 
many models: this has been verified either theoretically (e.g., for fBm or random wavelet series [33]), or 
numerically (e.g., for multifractal random walks or Levy processes, see [19]). The generic validity of the 
multifractal formalism has been proven in [46]. 

The p-leader multifractal formalism and its estimation performance are illustrated and studied in Sec¬ 
tion 5. 

3.3. p-leader Cumulants 

In [47] , the use of the cumulants of the log of p-leaders has been motivated by the fact that they enable 
to capture, with a small number of coefficients, most of the information contained in the spectra C^ p \h). 
Let Cm\j) denote the m-th order cumulant of the random variables log(£^ > ' ) ). Assuming that the moments 

(p) 

of order q of the p-leader t/f exist and 


via the Legendre transform 
(19) 


E 


o(p) 


— 2 -A p, (?)e 


(O 


where A 0 is the unit cube [0, l] d , one obtains that 

log (e[( 4 p) )«]) = log (e[(^V]) +C (p) (?)log(2 J 


( 21 ) 


7 











and, for q close to 0 


log (e [(#>)«]) =log (e^ 10 ®^]) = £ C&>(j)|^ 

m>l 


Comparison of (21) and (22) yields that the cumulants Cm\j) are necessarily of the form 

C%\j) = C%’V+c<£\og(V), 
and that p \q ) can be expanded around 0 as 

c (p) (,) = E^S- 

m> 1 


( 22 ) 


(23) 


(24) 


The concavity of (( p ' 1 implies that c% < 0. Note that, even if the moment of order q of is not finite, the 
cumulant of order m of log(£^) is likely to be finite and (23) above is also likely to hold. 

Relation (23) provides practitioners with a direct way to estimate the coefficients c'm in the polynomial 
expansion (24), termed the log-cumulants, by means of linear regressions of Cm\j) versus log(2 J ) [26]. 
Moreover, on condition that cij’ 1 ^ 0, the polynomial expansion (24) can be translated into an expansion of 


C^ p \h) around its maximum via the Legendre transform 

c^\h) = d+y % 


h — c 


( p )' 


l>2 


ml 


Jp) 


(25) 


with C 2 = cf\ C 3 = -Cg P \ C 4 = —cf^ + 3 ^ 55 -, etc., see [27]. 


3-4- Estimation and finite size effects 

In the computation of the multifractal formalism for data in applications (in contrast with the theoretical 
analysis of functions) one is confronted with the fact that only values of X sampled at finite resolution are 
available (in contrast with values for an interval in R d ). As a consequence, the infinite sum in the theoretical 
definition of p-leaders, cf. ( 8 ), is truncated at the finest available scale induced by the resolution of the data. 
Moreover, the liminf in the definition of the scaling function (18) cannot be evaluated. 

Our goal is here to compute the equivalents of (17), (18) and (23) for such finite resolution p-leaders and 
to use these expressions for identifying estimators for the p-leader scaling function (^ (q), the log-cumulants 
cfn and the p-spectrum C^ p \h). As a constructive model, we make use of a binomial deterministic wavelet 
cascade. For simplicity, we consider the univariate case d = 1. In Section 5, numerical results are provided 
that demonstrate that this simple deterministic model yields estimators with excellent performance for large 
classes of stochastic multifractal processes. 


3.4-1- p-leaders for deterministic wavelet cascades 

Let j < 0 denote the finest and j = 0 the coarsest available scales, respectively, i.e., j < j < 0. Let 
wo,wi > 0 and, by convention, let Cj— o,fc=i = 1. The coefficients ca of the deterministic wavelet cascades 
are constructed by the following iterative rule [9]. For j = 0, — 1,... ,j + 1, for each coefficient Cj^ at scale 
j two children coefficients are obtained at scale j — l by multiplication with ojq and ui 1 , Cj_i, 2 fc-i = <^o c j,k 
and Cj—i t 2 k — wiCj t k- At a given scale j < 0, there are hence 2 “i coefficients, taking the values 

Cj,k e {ujg n u}f :l ~ n , n = 0,..., -j}. 

The wavelet structure function (17) is therefore given by 

2 ~* 

s c (j, q ) = 2? J2( Cj , k y = 2 J (wg + = 

k =1 


Q\ -J 
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from which we identify the wavelet scaling function (7) of the cascade 


r](q) = l-log 2 (wg+u;?). 

Let p > 0 be such that rj(p) > 0. The restricted p- leaders, denoted fi p l , are defined by replacing 3A with A 

in (8) and the corresponding structure function, denoted S~ p \j,q-,j), by replacing i^ with in (17). It 
can be shown that structure functions with restricted p-leaders yield quantities equivalent to (17) so that 
the corresponding scaling functions (18) coincide, see [20]. Then 


m 





For infinite resolution j —> — oo, this expression boils down to = Cx ( t-a-iw ) ” ■ this case, the 

structure function reads 


?(p) 


00 ) = 2 3 ^{cj,k 

k=1 


1 


1 — 2 _, 7(p) 


_ 2^h(9) 


1 


1 _ 2~vip) 


SL 

v 


from which we identify the p-leader scaling function (18) of the cascade 

C (p) {q) = r)(q). 


However, for finite resolution j > — oo we have = c\ ^ ^ P and hence 


S?U, <,-,])= 2*°"^ j . 


(26) 


Similarly, if Cm' 1 (j) denotes the m-th order sample cumulant of log(f’^j) = log(cA) + ^ log ), 

then 


C[ p \j ) = c[ p ' 0) + c {p) log(2*) + - log 


l _ 2~U~j+ 1 )v(p) \ 


1 _ 2 ~v(p) 


) 


(27) 


while Cm\j ) = Cm\j) = Cm' 0 ' 1 +Cm^ log(2 J ) for to > 2 . Comparing (26) and (27) with (18) and (23) shows 
that the truncation of the definition of p-leaders (8) at finite scale induces an additional scale-dependent 
term, parametrized by j and p(p), in the p-leaders structure functions and first log-cumulant. 


3-4-2. Estimation 

Assuming that p(p) is known, the expressions (26) and (27) enable us to define estimators for C^ p \q) 

(p) 

and Cm as linear regressions in logarithmic coordinates 


3 2 


3 = 31 


&\q) = E w i ( lo §2 (s^\j,q)) - ^log 2 (l - 2 -C-1+1MP) 

J+i )v(p) \ 
1 _ 2 ~v(p) I 


(p) 1 v ' (r<r)r\ 

Cl = T —- 2^ Wj ( (j) — - log | 


log(2) 


m>2. 

V ' 3=31 


3 = 3 1 
h 


(28) 

(29) 

(30) 
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Estimators for the Legendre spectra can be defined in a similar way for the parametric development 
C{q), h{q) described in [48]. Note that the expressions (28) and (29) differ from those employed in the 
standard wavelet leader setting [19, 32, 26, 27, 21] due to the additional scale-dependent terms identified 
in (26) and (27). However, since rj(p) —* oo when p —> oo, these expressions are equivalent to the wavelet 
leader ones. In practice, the unknown function q[p) in (28) and (29) is replaced with the estimate 


J2 

v(p) = log 2 ( S cti’P)) ■ ( 31 ) 

3 = 3 1 

In the above expressions, j\ and j 2 are the finest and coarsest scales, respectively, over which the estimation 
is performed. The linear regression weights wj have to satisfy the constraints [T/j] jwj = 1 and Wj = 0 

and can be selected to reflect the confidence granted to each log 2 (iS 1 ^^(j, q)) (or Cm\j )), see e.g. [49]. In 
the numerical experiments reported in Section 5, following [49], we perform weighted linear fits; alternative 
choices have been reported in [26]. 

4. Multifractal Detrended Fluctuation Analysis 

4-1. Detrended Fluctuation Analysis for the estimation of the self-similarity parameter 

Detrended Fluctuation Analysis (DFA) was originally proposed for the estimation of the self-similarity 
parameter H for fractional Brownian motion (fBm) [50, 51]. Therefore, DFA has essentially been stated and 
studied in the univariate setting, d = 1. In essence, it relies on the observation that the Holder exponent 
of a sample path of fBm is a constant function h(t) = H. From the definition of the Holder exponent (cf., 
(12)) a natural multiresolution quantity T(a,t ) emerged 

Td(a,t) = \X(t) — Pt t a,N P (t)\ (32) 

where Pt, a ,N P is a polynomial of degree Np that is obtained by local fit to the data in a window of size a 
(note that, alternatively, the use of a moving average model for Pt ia ,N P has been proposed more recently, 
cf., e.g., [52]). From discrete versions T d (a, k) of these multiresolution quantities, structure functions 

S d {a, 2) = —Y j T d {a,k) 2 

are computed, and the parameter FI is classically estimated by linear regression of log 2 S d (a. 2) versus log 2 a. 
The performance of this procedure for the estimation of FI has been studied and compared to that of others, 
notably those based on wavelets, in various contributions (cf. e.g., [53, 54]). 

4-2. Multifractal extension 
4-2.1. Natural extension 

To extend DFA to multifractal analysis, a straightforward choice could have been to construct a multi¬ 
fractal formalism based on T d (a,t) (cf. [55] for a preliminary attempt): 

S d {a,q)=—YT d (a,ky ~C q a^,a-^ 0. 
n a 

a k 

However, it is now well understood that the Legendre transform of (, d (q) yields a poor upper bound for 
the multifractal spectrum V(h). One reason for this is that the values of T d concentrate around 0 and can 
thus not be raised to negative powers (note that this is also the case for increments or wavelet coefficients, 
leading to similarly poor estimates for T>(h)). From a theoretical point of view, this poor performance 
can be understood in the light of the fact that the definition of the Holder exponent (12) leads to the 
incorrect intuition that it must be tied to increments, X(t + a) — X(t), or “improved” increments, T dl as 
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multiresolution quantities. However, recent theoretical contributions (cf., e.g., [19, 32, 26, 56, 21]) show that 
correct multiresolution quantities associated with the Holder exponent are the oscillations 

0{a,x) = sup \X(u) — X(v)|. 

u,v£B(x,a) 


A relevant multifractal formalism can be constructed for 0{a,x) as long as the smooth parts of the data 
have Holder exponents strictly smaller than 1. Otherwise, higher order oscillations must be used. Note that 
in a wavelet framework, wavelet leaders play for wavelet coefficients the same role as oscillations play for 
increments. 

4-2.2. L 2 -norm formulation 

The use of oscillations 0{a,x) instead of Td would thus have been a natural way to extend DFA to 
multifractal analysis. However, the multifractal extension proposed by Kantelhardt et al. followed a different 
path: In the seminal contribution [23], they proposed the following original multiresolution quantity 

T m fd{a, k) = ^ ^ \X(ak + i) - Pk,a,N P (i)\ 2 ^j ,k = !,•••>“ (33) 

and constructed a multifractal formalism with T m fd(a,k), the so-called multifractal detrended fluctuation 
analysis (MFDFA). Here n denotes the number of available samples and Pk,a,N P is the same polynomial as 
in (32). Defining structure functions 

n/a n/a / a \ § 

Smfd(a,q) = -YT mfd (a,k) q = - Y ( - Y \X(ak + i) - Pk, a ,N P {i )| 2 ) 
n k n k\ a ti ) 

and the corresponding scaling function 


r , \ y ■ ( iog (S m fd{a,q)) 

C mfd(q ) = lim inf -pyr- 

j->-oo log(a) 

yields, via a Legendre transform, the MFDFA multifractal formalism 

£mfd{h) = inf (d + qh - Cmfd{q)) > T>{h). 

q€ R 

It was numerically compared against the wavelet leader and WTMM multifractal formalisms for multiplica¬ 
tive cascade-type multifractal processes, see [23] for the original contribution. Despite the fact that there has 
been no theoretical motivation for the use of (33), the MFDFA formalism was found to perform very satis¬ 
factorily and is commonly used in applications (cf., e.g., [57, 58, 59, 60, 61, 62]), thus empirically justifying 
a posteriori the choice of T m f d as a relevant multiresolution quantity for multifractal analysis. 

4-3. p-exponents, p-leaders and MFDFA 

Comparing the definition of the multiresolution quantity T ^ in (11) for the p-exponent to that of T m f d 
in (33) for MFDFA clearly shows that the latter mimics the former, for p = 2, with a discretized setting of 
the continuous integral in (11), and with the theoretical Taylor polynomial P Xo replaced by a data-driven 
locally adjusted polynomial Pk,a,N P - Therefore, MFDFA can a posteriori be interpreted, and theoretically 
grounded, in the framework of p-exponent analysis: While MFDFA was originally associated to the analysis 
of local regularity as measured by the Holder exponent, the present contribution thus shows that T m f d 
must be related to a p-exponent based characterization of local regularity with p = 2 and not to the Holder 
exponent. 

The following paragraphs provide a detailed theoretical and practical comparison of the MFDFA and 
p-leader frameworks. 
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Choosing p. In the MFDFA method, p is arbitrarily set to 2, while p-exponents in (11) are theoretically 
defined for p £ (0,po). There are practical benefits stemming from varying p, as summarized in Section 2.3 
and discussed in detail in the companion paper [35, Section 4]: Notably, the analysis requires to choose p 
such that rj(p) > 0 ; moreover, the use of various values p < Po may help practitioners to reveal the fine 
local singularity structure of data (see Section 5.3 for a numerical illustration). 

Local regularity and integration. In numerous real-world data, notably in biomedical applications 
(Heart Rate variability, fMRI), it is observed that h min < 0 (cf., e.g., [28]). This explains why MFDFA 
procedures, as detailed, e.g., in [23, 63, 64], always perform an integration of the data, X(t) -£ f* X(s)ds, 
as a preliminary step (hence implicitly assuming that for the integrated data, h min > 0). However, as long 
as r/(p) > 0 and h > —1/p, this preprocessing of data is not needed and may even constitute a drawback if 
data contain oscillating singularities, cf. [35, Section 4]. 

Time domain versus wavelet domain. The MFDFA method relies on a time domain implementation, 
i.e., T m fd is computed directly from the sampled time series X, while the p-leader formalism relies on 
wavelet projections. This has fundamental theoretical and practical implications which are detailed in the 
next paragraphs. 

Local regularity, Taylor polynomial and vanishing moments. A key property of the wavelet 
characterization is that it does not require the knowledge of the Taylor polynomial (cf. [21]) whereas 
MFDFA require some (heuristic) estimation of this polynomial (cf. a contrario [24, 25], where a variation 
of MFDFA that attempts to avoid the estimation of the polynomial is devised). Despite its being inspired 
from the definition of the Holder exponent, the polynomial Pt, a ,N P {t) is in practice obtained as the best fit 
of X(u) for u £ [t — a/2,t + a/2] for a priori chosen and fixed degree Np. The polynomial Pt } a,N P must be 
computed for each time position t and analysis scale a. It thus depends on a, while it should not depend on 
a in theory, cf., (11). Moreover, the order of the polynomial in (11) can depend on the time position t (cf., 
[21]) while the order Np in the MFDFA method is fixed. 

Another key difference in the use of p-leaders and T m fd lies in the fact that the former requires the 
choice of the number of vanishing moments N,p of the mother wavelet ipo, while the latter implies the choice of 
the degree Np of the polynomial Pt, a ,N P ■ Often, the parameters and Np are regarded as equivalent, yet 
this is only partially correct since the choices of N,p and Np are framed by different theoretical constraints: 
In order to recover the local power-law behavior T^\a,t) ~ Ca h ,a 0 for an isolated singularity with 
regularity exponent h , N,p must be larger than h (and can be set globally to be larger than the largest 
regularity exponent of X, cf., [35, Section 3.1], [19]), while it is required that Np < h by definition of the 
Holder exponent. 

Robustness to trends and finite size effects. Practically, the choice of and Np is often thought 
of in terms of robustness to trends (cf. e.g., [65, 50, 53, 66, 67, 68, 69]). 

Assuming that the data to analyze are actually corrupted by deterministic trends Z seen as noise, 
X + Z, it has been documented that increasing the number of vanishing moments of the mother wavelet Ay, 
diminishes the impact of the additive trends to the estimation of the scaling exponent [65, 53]: The possibility 
of varying N,p brings theoretical and practical robustness to analysis and estimation. The practical price 
paid for increasing consists of a larger number of wavelet coefficients being polluted by border effects 
(finite size effect) and can also decrease the coarsest scale at which data can be analyzed. This price turns 
out to be very low since coarse scales contain only few coefficients and hence have very little impact on 
estimated scaling exponents. Note that the finest scale that is available using p-leaders is independent of 
N,p and directly and only determined by the resolution of the data (i.e., if the data is sampled at At, then 
the the finest scale contains coefficients at rate 2At). 

In the MFDFA method, increasing Np also brings some form of robustness to additive noise, however, 
this far more depends on the nature of the trend [66, 67, 68, 69] (see also the numerical illustrations in 
Section 5.2.3). Yet, increasing Np has a more drastic practical consequence: The finest scales of data 
cannot be used practically since the best fit of X(u) for u £ [t — a/2, t + a/2] cannot be achieved until the 
number of samples actually available in \t — a/2, t + a/2] is substantially larger than Np: For a polynomial 
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of order Np, one needs at least Np + 2 data points and consequently, A j = [log 2 (lVp + 2)] — 1 fine scales 
that can be analyzed with p-leaders are not accessible for MFDFA. Losing fine scales is problematic from a 
statistical estimation point of view (since fine scales contain many coefficients) as well as from a multifractal 
analysis point of view (since it theoretically requires to estimate the scaling exponents in the limit of fine 
scales a —> 0). 

In summary, the choice of Np, is theoretically better grounded and practically much easier and less critical 
than that of Np. The polynomial subtraction entering MFDFA is thus a more intricate issue than it may 
seem, with little theoretical guidelines. 

Extension to higher dimension. Theoretically, MFDFA can be extended straightforwardly to higher 
dimensions, d > 1. In practice, the computation of local best fit polynomials becomes a real issue in higher 
dimension, and there are few attempts to extend MFDFA to dimension d = 2 [70, 71, 72, 73]. Along 
the same line, the use of polynomials of degree 2 or higher is uneasy (higher order polynomial fitting and 
multivariate integral numerical evaluation). In contrast, the p-leader formalism extends without difficulties 
to higher dimensions given that higher dimensional discrete wavelet transforms are readily obtained by tensor 
product of ID-mother wavelets (cf., e.g., [74]). Also note that the MFDFA method has larger computational 
complexity (of order O(nlogn), where n is the total number of samples) than the p-leader formalism (of 
order 0(n)), which quickly becomes an issue in higher dimension. 

Conclusions. This connection of MFDFA to p-leaders provides a theoretical grounding for the choice 
of Tm/d, which otherwise appears as a relevant, yet ad-hoc, intuition used to construct a multifractal 
formalism. It also makes explicit that MFDFA measures local regularity via the 2-exponent and not the 
Holder exponent. The p-leader formalism can thus be read as an extension (different p) and wavelet- 
based reformulation of MFDFA. The p-leader and MFDFA formalisms are compared in terms of estimation 
performance and robustness to trends in Section 5.2.2 below. 

It is also of interest to note that p-exponents and p-leaders, on one side, and MFDFA, on other side, 
were proposed independently and in different fields: In the Mathematics literature for the former around 
year 2005 [34, 40], in the Physics literature for the latter around year 2002 [23]. This is, to the best of our 
knowledge, the first time that these two notions are connected, related and compared. 

5. Illustrations and estimation performance 

We numerically illustrate and validate the proposed p-leaders multifractal formalism and compare it 
against the leader and MFDFA formalisms. To this end, the formalisms are applied to independent realiza¬ 
tions of synthetic random processes with prescribed multifractal p-spectra, and their respective estimation 
performances are studied in detail and compared. For the sake of exhausitivity in comparisons, processes 
whose p-spectra do not depend on p as well as processes whose spectra vary with p are used. Also, situations 
in which p-exponents all collapse with the Holder exponents are considered as well as examples where this is 
not the case. The p-leader multifractal formalism is furthermore numerically validated in higher dimensions 
by application to synthetic (2D) images with prescribed multifractal p-spectra. 

The WTMM formalism has already been compared independently against the wavelet-leader formalism 
[33, 32] and MFDFA [23, 75] and is thus not re-included in this study. 

Results illustrate the benefit of the extra flexibility of varying p in the p-leader multifractal formalism, 
both in terms of estimation performance and for evidencing data whose multifractal p-spectra are not p- 
invariant. 

Sample paths of all processes were numerically simulated by Matlab routines implemented by ourselves, 
available upon request. 
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5.1. Random processes with prescribed multifractal p-spectra 

5.1.1. Fractionally differentiated Multifractal Random Walk 

The multifractal random walk (MRW) has been introduced in [76] and is defined as 

n 

X{k) = ^2 G H (k) exp (w{k)) 

k=1 

where G#(fc) are increments of fractional Brownian motion with parameter H > 1/2 and u is a Gaussian 
process that is independent of Gh and has covariance cov(w(fci), ui(k 2 )) = Ain ^ \k 1 - C k 2 \+'^ ) when \ki — 
k '2 < L , and 0 otherwise, with A > 0. By construction, MRW has stationary increments and mimics the 
multifractal properties of Mandelbrot’s multiplicative cascades [9]. It has been chosen here as an easy-to- 
simulate member of this widely used class of multifractal processes. 

By means of fractional integration of negative order s < 0 of X (in practice fractional differentiation of 
positive order v = —s, cf., [35, Definition 3]), we obtain sample paths X^ v > with different values for the 
critical Lebesgue index po- MRW contains only canonical singularities (cf. [35, Definition 4]), hence fractional 
differentiation results in a pure shift of their multifractal p-spectra by v to smaller values of h. Moreover, 
since canonical singularities are p-invariant (cf., [35, Theorem 2]), the multifractal p-spectra of collapse 
for all p < po and are given by 


V ( f>\h)=V u {h) = 


1 + £ 2 . 
X ' 2 


(^j 


for h p £ 
otherwise 


ci,*, — \/~ 2C2, Cl,*, + y/— 2C2 


(34) 


where ci „ = H + A 2 /2 — v, c -2 = —A 2 . Furthermore, c m = 0 for all m > 3. The wavelet scaling function is 
given by 

, , _ j (H + A 2 /2 - u)p + ffp 2 for 0 < p < yJ—2 /c 2 

|(i7 + A 2 /2 - v)y/-2/c 2 - 1 - y/-2c 2 + c 2 p for p > ^/—2/c 2 . 

By elementary calculations, evaluation of condition (13) yields 


oo for v £ [0, H + A(| — \/2)] 

p 0 = { 1/ (v — H — A(| — \/2)) for v £ [H + A(^ — \/2), H + A(f - ^)] 

2(H + - v)/\ 2 for v £ (H + A(| - ^), H + %] 


(35) 


and p(p) < 0 for any p>0ifv>H + A 2 /2. 


5.1.2. Lacunary wavelet series 

Lacunary wavelet series (LWS) X a depend on a lacunarity parameter p £ (0,1) and a regularity 
parameter a £ R. At each scale j < 0, the process has a fraction of exactly 2 _TU nonzero wavelet coefficients 
on each interval [1,1 + 1) (l £ Z), at uniformly distributed random locations, and whose common amplitude 
is This process was introduced in [77] when a > 0 (i.e., suited to Holder exponent based multifractal 
analysis) and in [78] for the general case a £ M (thus requiring the use of p-exponents). The use of LWS 
here is motivated by the fact that they contain singularities that are not p-invariant: Indeed, it is shown in 
[78] that almost every point is a lacunary singularity (cf., [35, Definition 7]). For a > 0, p £ (0,1), LWS are 
bounded and their p-spectra Vp > 0 are given by [78] 


V (p \h) 


H+l/p 
^ a + 1/p 


h £ 



—oo otherwise. 


(36) 
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5.1.3. Simulation setup 

We generate Nmc — 500 independent realizations of N = 2 16 samples each of MRW (with param¬ 
eters H = 0.72, A = \/0.08 and fractional differentiation v £ {0, 0.4, 0.6, 0.7, 0.73}, yielding p 0 = 
{+oo, 25, 4, 1.5, 0.75}, cf., (35)) and LWS (with parameters a £ {0.2, 0.3}, p £ {0.7,0.8}, po = +oo). For 
each realization, we compute the Legendre spectra C^ p \h p ) in (19) and C m fd{h ) as well as the log-cumulants 
Crn for m = 1, 2, 3 (cf., (23-24)). We adhere to the convention that the finest available dyadic scale is labelled 
by j = 1. Estimates are computed using (28-30) for dyadic scales from j\ = 4 to the coarsest available scale 
j*2 (j -2 = 13 for p-leaders due to border effects and j 2 = 15 for MFDFA) with weighted linear regressions 
(bj = nj , see [49]). For p-leaders, a Daubechies’ wavelet with TV,/, = 2 vanishing moments is used, and 
consistently the degree of the polynomial for MFDFA is set to Np = 1. Furthermore, the p-leader estimates 
are calculated for the set of values p £ {j, 1, 2, 4, 5, 8, 10, +oo} (where p = +oo corresponds to 

wavelet leaders). 

5.2. p-invariant p-spectra and negative regularity 

We use fractionally differentiated MRW here because it enables us to study negative regularity and to 
compare the respective performance for different values of p since its multifractal p-spectra are p-invariant. 

5.2.1. Numerical illustrations of multifractal p-spectra 

We illustrate and qualitatively compare the p-leader, leader and MFDFA multifractal formalisms, based 
on their Legendre spectra C^ p \h) and C m fd(h). Averages over independent realizations are plotted in Fig. 
1 (colored solid lines with symbols), together with the theoretical spectra (34) (black solid lines and shaded 
area) and the respective theoretical bounds (16) for (colored dashed lines). 

Sample paths. Fig. 1 (left column) plots representative examples of sample paths with, from top to 
bottom, increasing value of v (and, hence, decreasing regularity and po). Visual inspection of the sample 
paths indicates the practical benefit of the use of model processes with (potentially negative) regularity: 
While (positive only) Holder regularity implies relatively smooth sample paths, the use of (negative) p- 
exponents provides practitioners with a rich set of models for applications with highly irregular sample 
paths, with a continuously rougher appearance as p-exponents take on smaller and smaller values (cf. first 
column of Fig. 1). 

2-leaders and MFDFA. In Fig. 1 (center column) Legendre spectra obtained with MFDFA and p-leaders 
with p = 2 are plotted together with the theoretical 2-spectra 

- The estimates C^ 2 \h) and C m fd{h) are observed to be qualitatively equivalent when po is large. 
However, for negative regularity (2 < p 0 <C +oo), the Legendre spectra C m fd(h) obtained with MFDFA 
only partially capture the theoretical spectra for negative values of h and appear shifted to larger values 
of h with respect to the theoretical spectrum V l ' 2 \h). 
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Figure 1: Fractionally differentiated MRW with p 0 = {+oo, 25, 4, 1.5, 0.75} (from top to bottom 
row); the top row plots MRW without fractional differentiation: Single realizations (left column), theoretical 
spectra T>(h) (black solid line and shaded area), estimates C m fd{h ) and C (2 ^ ( h ) (center column) and estimates 
(right column). The dashed lines indicate the theoretical bound (16). 
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- In contrast to MFDFA, the 2-leaders formalism clearly provides excellent estimates for ( h ) for any 

of the processes for which po > 2. 

We conjecture that the bias (shift towards positive values of h) of C m fd{h) for processes with negative 
regularity is caused by finite size effects similar to those analyzed for p-leaders in Section 3.4. 

p-leaders for p ^ 2. In the right column of Fig. 1, averages of p-leader estimates C^ p \h) are plotted for 
several values of p and compared to the theoretical p-spectra V^ih) and the theoretical bounds (16). 

- Clearly, the p-leader multifractal formalism provides excellent estimates for the theoretical spectra 
T>l P \h) for p < po, hence validating the proposed formalism. 

- The estimates are of excellent quality also for p < 1; notably, the choice p = 1/2 < 1 enables to 
correctly recover the p-spectrum T>l P \h) for with po = 0.75 (Fig. 1, bottom row), which is not 
possible for p > 1 (and, hence, neither with the MFDFA method). 

- When p > po, the estimates C^\h) are tangent to the theoretical bounds (16). Consequently, they 
are shifted to larger values of h. with respect to the theoretical spectrum v[f\h) and hence biased. 
This is visually most striking for the case p = +oo (i.e. for classical leaders associated with Holder 
exponents), for which estimated spectra are constrained to positive values of h. This phenomenon will 
be investigated in a forthcoming study (see [79] for preliminary results). 

- Finally, in consistency with [35, Theorem 2], the spectra C( p \h) coincide for all p < p 0 for fractionally 
differentiated MRW. 


5.2.2. Estimation performance 

We proceed with a quantitative analysis of the estimation performance of the p-leader and MFDFA 
multifractal formalisms, respectively. To this end, we assess the estimation performance for the log-cumulants 
cfn for to = 1, 2, 3,4 (cf., (23-24)) based on their root mean squared error (rmse), defined as 



where ( • )n mc stands for the average over independent realizations. 

Results for fractionally differentiated MRW with v e {0, 0.4, 0.6, 0.7} (po = {+oo, 25, 4, 1.5}) 
are plotted in Fig. 2 (top to bottom rows, respectively); the logarithm (log 10 ) of rmse values of p-leaders 
are plotted as a function of p (solid red lines with squares), those of MFDFA are plotted with a red circle. 
Distributions of the estimates, after subtraction of theoretical value, are shown in black boxplots for p-leaders 
and blue boxplots for MFDFA. The values for p on the right of the vertical red dashed lines are larger than 
Po, p > Po- Since c,? = c m , we omit the superscript below. 

Estimation of c±. Rmse values and estimates for the first log-cumulant are reported in the left column of 
Fig. 2 and yield the following conclusions. 

- For p < po, the p-leader multifractal formalism systematically yields better estimation performance for 
small values of p than for large values of p; the improved performance for small values of p is induced 
by reduced standard deviations of estimates, resulting in considerable rmse gains of up to a factor 2 
for p = 1/2 over the value of p which is picked closest to po- 

- For p > po, there is a sharp increase in rmse due to a systematic bias. Indeed, ci captures the position 
of the maximum of the p-spectra, which is shifted to larger values of h as compared to the theoretical 
spectrum when p > po, as discussed in Section 5.2.1. 

Estimation of c m for m > 2. The second, third and fourth column of Fig. 2 plot rmse values and 
estimates for C2, C3 and C4, respectively, and yield the following conclusions: 


17 



p 0 = +oo, h mm = +0.36 
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Figure 2: log-cumulants of fractionally differentiated MRW with p 0 = {+oo, 25, 4, 1.5} (from top 
to bottom). Red: log 10 of rmse of p-leaders (solid lines and squares), as a function of p, and for MFDFA 
(full circles), for c m , m = {1, 2, 3} (left to right column). Black and blue: boxplots of the error (c m — c m ) 
for p-leaders (black) and MFDFA (blue). The left and right y-scales correspond, respectively, to rmse and 
boxplots. Rmse is shown in the same y-scale across the rows. 

- A systematic benefit in choosing small values for p in the p-leader multifractal formalism is observed. 
Indeed, rmse values decrease sharply and pronouncedly for p < 4, and the smallest rmse values are 
obtained for small values for p (p = 1/2 or p = 1). 

- Unlike for ci, choosing p > po does not significantly alter estimation performance for c m for m > 2. 
The interested reader is referred to [79] for further details. 

Comparison with MFDFA. 

- The MFDFA and 2-leader formalisms have similar performance for the estimation of C\ as long as 
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Figure 3: MRW with non-polynomial trend ( H = 0.72, A = a/0.08, u = 0), for different detrending 
powers: = 3, Np = 2 (top row), TV/, = 4, Np = 3 (center row), TV/, = 5, Np = 4 (bottom row). Red: 

rmse for p-leaders (solid squares), as a function of p, and for MFDFA (full circles). Black and blue: boxplots 
of the error (c m — c m ) for p-leaders (black) and MFDFA (blue). The left and right y-scales correspond, 
respectively, to rmse and boxplots. Rmse is shown in the same Y-scale across the rows. 


Po = +oo. As soon as the support of the p-spectrum includes negative values for h and po < oo, 
MFDFA produces estimates for ci that are biased. This is consistent with the observations of Section 
5.2.1 where the C m fd{h ) are found to be shifted to larger values of h. 

- For the estimation of c m for m > 2, the MFDFA formalism yields similar performance as the p-leader 
formalism for moderately small values p sa 2. 

These results lead to the conclusion that, for data containing p-invariant singularities only, it is beneficial 
to choose a small value of p in the analysis. Note that the p-leader formalism with moderately small values 
for p, e.g. p < 4, significantly outperforms the current state-of-the-art wavelet leader formalism (p = +oo) 
which yields up to 50 percent larger rmse values. The MFDFA and p-leaders formalisms have comparable 
performance for the estimation of c m for m > 2 and also for C\ as long as po = +oo. Yet, MFDFA estimates 
of Ci are biased when the data are characterized by negative regularity exponents. 

5.2.3. Robustness against non-polynomial trend 

Here, we study the robustness of MFDFA and the p-leader formalism for estimation of the log-cumulants 
c m for m = 1, 2,3,4 when a G°° but non-polynomial trend of the form 

T (t) = 100(f+ 1/100)“ 1/2 , ie[0,l] 

is added to MRW with parameters specified as in Section 5.1.3 and v = 0. We use Daubechies’ wavelet with 
Ni[, = {3,4,5} vanishing moments for p-leaders and, correspondingly, polynomials of degree N P = {2,3,4} 
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a=0.3, r|=0.8 


a=0.2, »)=0.8 


Figure 4: LWS for (a,rj) = {(0.2, 0.7), (0.2, 0.8), (0.3, 0.7), (0.3,0.8)} (from top to bottom row): Single 
realizations (left column), theoretical spectra V(h) (right column, dashed lines) and estimated spectra C(h ) 
(right column, solid lines). 


for MFDFA. Note that with the largest choice Np = 4 for the degree of polynomial, the finest available scale 
for MFDFA is j = 4; unlike p-leaders, finer scales cannot be used for estimation with MFDFA, cf., Section 
4.3. In order not to penalize MFDFA in the comparisons, we perform linear regressions from scale j\ = 4 to 
the coarsest available scales (j? = 13 for p-leaders due to border effects and j -2 = 15 for MFDFA). 

Rmse values and estimates for c m , m = 1,...,4, are plotted in Fig. 3 for = {3,4,5} (top to 
bottom rows, respectively). Clearly, the estimation performance of MFDFA is severely degraded by the 
non-polynomial trend: Even with the polynomial of highest degree considered here (Np = 4), rmse values 
for C2, C3 and C4 for MFDFA are up to one order of magnitude larger than in the absence of the trend (cf. 
Fig. 2, top row). 

In contrast, the rmse values for p-leaders with N^. — 4 are very close to those obtained in the absence 
of the trend. This indicates that the wavelet transform underlying p-leaders is considerably more effective 
in removing the impact of the trend on the higher order statistics of the multiresolution quantities than the 
empirical polynomial-fitting procedure employed by the MFDFA method. 
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5.3. p-spectra and lacunary singularities 

We illustrate the p-leader multifractal formalism for the estimation of the multifractal p-spectra of LWS, 
which contain lacunary singularities at almost every point and are hence not p-invariant. Averages of 
£( p )(/i) over independent realizations are plotted in Fig. 4 (colored solid lines with symbols), together 
with the theoretical spectra (36) (colored dashed lines) for four combinations of the parameters (a,rj) with 
a £ {0.2, 0.3} and rj £ {0.7, 0.8}. Since the MFDFA method cannot reveal the difference in the p-spectra 
because it is limited to p = 2, it is omitted in Fig. 4. 

As expected (cf., (36)), the numerical estimates of the p-spectra are not invariant with p but reproduce the 
evolution with p of the theoretical spectra T>( p \h): The larger p, the more the upper limit of the support of 
the spectra are shifted towards the point a. The positions of the mode of C^ p \h) slightly underestimate those 
of the true spectra, and the Hausdorff dimension of the leftmost point of the spectra are poorly estimated. 
Yet, the C^ p \h) qualitatively reproduce the theoretical spectra satisfactorily well and, in particular, clearly 
and unambiguously reveal the lacunary nature of the sample paths. 


5.4- Images: Canonical Mandelbrot Cascades 

As mentioned above, MFDFA has barely been used for images (except for the attempts in [70, 71, 72, 73]), 
while the p-leader multifractal formalism extends in a straightforward manner to higher dimensions. Here, 
we illustrate this point and apply the p-leader multifractal formalism to synthetic multifractal images. 
Canonical Mandelbrot Cascades (CMC). The construction of multiplicative cascades of Mandelbrot 
(CMC) [9] is based on an iterative split-and-multiply procedure on an interval; we use a 2D binary cascade 
for two different multipliers: First, log-normal multipliers (CMC-LN), W = 2~ u with U ~ Af(m, 2m/ ln(2)) 
a Gaussian random variable; Second, log-Poisson multipliers (CMC-LP), W = 2 7 exp (1ii(/3)7Ta), where n\ is 
a Poisson random variable with parameter A = — ■ We use fractional integration of order a = 0.2. CMC 

contain only canonical singularities, hence fractional integration results in a pure shift of their multifractal p- 
spectra by a. Their multifractal p-spectra hence all collapse due to the p-invariance of canonical singularities. 
For CMC-LN, the multifractal p-spectrum is given by 

D^ p \h) = D(h) = 2 — --—-—, with C\ = m + a, C2 = —2m, c m = 0 for all m > 3. 

4m 

The expression for the multifractal p-spectrum of CMC-LP is 


V(h) =2 + 


- 1 


—a + 7 + h 


(—a + 7 + /i )(/3 — 1) _ 

7 In (3 


with ci = a + 7 — lj and all higher-order log-cumulants are non-zero: c m = (— ln(/3)) m . We set 

m = 0.04, p = 0.8395 and 7 = 0.4195, yielding ci = 0.24 and c-i = —0.08 for both cascades, and C3 = 0.014 
for CMC-LP. 

Illustration of multifractal p-spectra. Averages over 100 realizations of p-leader Legendre spectra 
C^(h) 2D CMC of side length N = 2 11 are reported in Fig. 5 (bottom) for CMC-LN (left) and CMC-LP 
(right), single realizations of the random fields are plotted in Fig. 5 (top); we use the scaling range j 1 = 3 and 
j -2 = 8 and tensor-product Daubechies’ wavelet with TVy, = 2 (cf. [74]). The observations and conclusions are 
similar to those obtained in the previous subsection for (ID) signals (MRW); indeed, the p-leader estimates 
C^ p \h) (solid lines in color, symbols) enable to correctly recover the theoretical p-spectrum D^ p \h) as soon 
as the condition p < p 0 is fulfilled. 


6. Discussions, Conclusions and Perpectives 

The present contribution has developed and analyzed a novel multifractal formalism based on new local 
regularity exponents, the p-exponents, and their corresponding multiresolution quantities, the p-leaders, 
that have been theoretically defined and studied in the companion paper [35]. This new formulation of 
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Figure 5: 2D Mandelbrot cascade. Log-Normal (left) and Log-Poisson (right) multipliers. Top: Single 
realization. Bottom: V(h) (black solid line and shaded area), C^ p \h) for p = +oo (black, cross), p = 1/2 
(blue, circle), p = 2 (magenta, square), p = 4 (red, triangle), p = 16 (green, diamond) and theoretical limits 
for multifractal p-spectra (dashed). 


multifractal analysis generalizes the traditional Holder-exponent-based formulation in several ways. First and 
foremost, it naturally allows to perform the multifractal analysis of functions that have negative regularity 
or, equivalently, that are not locally bounded (but instead belong locally to L p ). This allows to avoid the 
commonplace a priori (fractional) integration of the data and the pitfalls it entails. 

Moreover, the dependence on the parameter p provides an important and rich information concerning 
the characteristics of the singular behavior of data. When the multifractal spectra differ for different valid 
values of p, this clearly indicates the presence of complex singular behaviors on the data, known as lacunary 
singularities. This information is a distinctive feature of p-exponents and is not accessible with previous 
tools. 

This contribution has also made a clear connection between the p-leader multifractal formalism and 
MFDFA, a related technique that is widely used in applications. This has allowed to provide a theoretical 
grounding to MFDFA, which remained a useful but ad hoc intuition. Also, it has brought to light that 
MFDFA actually measures the 2-exponent, rather than the Holder exponent as had been assumed previously. 

Numerical simulations on synthetic multifractal processes have shown that the p-leader formalism, for 
small values of p, benefits from significant improvements in estimation performance, compared to wavelet 
leaders and MFDFA. It is important to emphasize that, even though p-leaders were originally introduced 
with the goal of analyzing negative regularity, they show better estimation performance than wavelet leaders 
even for processes that have positive regularity only. 

To conclude, the benefits of the p-leader multifractal formalism over existing formulations are threefold: 
First, it allows the analysis of negative regularity; second, it shows better estimation performance; and third, 
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the dependence of the estimates on the parameter p provides richer and more detailed information on the 
characteristics of the singularities that can be found in the data. This suggests that the p-leader multifractal 
formalism, for small values of p , should be preferred for practical multifractal analysis. 

A Matlab toolbox implementing estimation procedures for performing p-multifractal analysis will be 
made publicly available at the time of publication. 
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Appendix A. Proofs 


Proof of Proposition 1: We start by considering the case p > 1, in which case, we will prove the 
slightly sharper result that (16) holds as soon as / £ L p . Let H > —d/p be given. Let Eh denote the set of 
points where the p-exponent of A' is smaller than H. If Xg £ Eh, then there exists a sequence r n —> 0 such 


that 


I B(x 0 ; 


\ 1 /P 

) | f(x)\ p dxj > r”, 


so that there exists a sequence j n —> — oo such that the dyadic cubes Xj ri (xg) satisfy 


'3\j n (x 0 ) 


\f{x)\ p dx > C2^ d+Hp ) 


(pick the smallest dyadic cube Xj n (x 0 ) such that B(x 0 ,r n ) C 3Aj ii (x 0 )). One of the 3 d cubes of width 2 P 
that constitute 2>Xj n (which we denote by pj n ) satisfies 


' Pin (* 0 ) 


\f{x)\ p dx > C3~ d 2^ d+Hp '> 


We consider now the maximal such dyadic cubes, of width less than a fixed e, satisfying this inequality for 
all possible xq, and we denote by 1Z this collection. Then, since maximal dyadic cubes necessarily are 2 by 
2 disjoint, 

c 2^ d+Hp) I dx - c 

fj,eiz p 

where p is of width 2 J . If a; £ Eh, then x belongs to one of the 3/x, and therefore to the ball of same center 
and radius 3d2# Since r n > C 2 J ’*, we have obtained a covering of Eh by balls of radius at most 3 de such 
that 

^ 2diam(B(x,r)) d+Hp < C, 

and the result follows for p > 1. 

We now consider the case the case p < 1. Hypothesis p(p) > 0 means that 


3C, e > 0 : Vj < 0, 2 dj ^ \cf k \ p < C2 £pj . 

i,k 


(A.l) 


If h p ( xq) < H, then there exists an infinite sequence of dyadic cubes A which contain xg and such that 


H l c v T 2" d(i_j,) > 2 Hpj . 

j'<j, A'C3A i =1 
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We now consider the maximal cubes of width less that e and which satisfy this condition. This yields a 
covering of Eh by dyadic cubes A which are 2 by 2 disjoint. Denote by Aj the cubes of this covering which 
are of width 2 J , and by Nj the cardinality of Aj. On one hand 

E E \°y\ P 2 _d0_/) > Nj2 Hpj 

A j'<j, A'C3A i=l 


(where the sum over A is taken on all dyadic cubes of width 2 J "). On the other hand, the left side is equal to 


2 d —1 


3d E E E l^l p 2-^ 

j'<j \k’e1 d i —1 

But (A.l) implies that the term between parentheses is bounded by 2 — 2 epj ; thus we obtain that 

Nj2 Hpj < C2~ dj , 

which implies that dim(.Eff) < d + Hp, and the result follows for p < 1. 
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